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The processes of Coulomb gas ordering in 3D layered system are studied by 
means of Brownian dynamics approach. It is found that at different densities 
of the carriers the 3D lattice of charges as well as new specific structures 
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are possible in the system. In particular, at some small density the particles 



inside the layers can associate into droplets that collectively repulse between 
r**"* ' neighbouring layers. These droplets possess local stripe structure which orders 



spontaneously along arbitrary chosen direction. At higher densities a specific 



Q\ ' ordering of the charges into the tetragonal-like or hexagonal-like structures is 

On' 

observed visually and described numerically. Specific "pairing" of the charges 

from different layers plays an essential role in formation of all above structures. 

'"O 

It is essentially many-body effect in 3D system which can be a reason for 

o" 

^ ' unusual properties of layered crystals. 
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I. INTRODUCTION 

The theory of doped Mott insulator systems have shown variety of electronic ground 
states from electronic liquid crystal phases!! to the usual metallic state that depend on dop- 
ing level. The complex character of phase diagram is caused by the long-range Coulomb 
interaction between electrons placed in an antiferromagnetic background. Besides, the well 
defined layered structure is a common feature of these systems at all variety of other prop- 
erties. Both these features are essential difficulties that obstacle to construct correct theory 
here up to the moment. 

Therefore, the majority of the works in the field are the studies of the two-dimensional 
systems with quite complicated interactions. The separation of the charge and spin variables 
is an additional simplification hereall Sometimes such models are not quite realistic and can 
not be applied to describe the systems without magnetic ions. For example, some systems 
without magnetic ions, like Bai^rC^-BiOa, have shown a remnant of charge ordering above 
T c -temperature of superconductive phase transition for the optimally doping value of x=0.4. 

In present article we show that the density stratification (or phase separation) accompa- 
nied by a collection of classic phases (like Wigner crystal or stripe-ordering) appears already 
in a very simple model of three dimensional layered crystal with long-range Coulomb inter- 
action among spinless charges. 

One deals here with a case that is extremely difficult for the theory. From one side the 
well pronounced anisotropy and layered structure do not give a chance to limit ourselves by 
simple mean-field approximation which is suitable for analytic theory. From another one, 
the essentially many-body and 3D character of the system makes it difficult to simulate its 
behavior directly. 

However, the development of the computer technology continuously lowers the extent 
of the simplification that is necessary for numerical simulation. As a rule, the simulation 
in two or three dimensions is performed separately. Keeping in mind an application to 
above layered structures, it is interesting to simulate just intermediate case. Let us suppose 



that the equally charged spinless particles move and interact in 3D space with a positive 
background that includes a periodic potential along one of the dimensions which we choose 
as z-direction. 

First of all, let us discuss briefly qualitative effect of "layer structure" . At some conditions 
the particles should be found as localized completely in a vicinity of the equidistant layers. 
It generates a strong impact from the 2D behavior into collective behavior. The particles 
tend to produce a hexagonal (Wigner) crystal in each layer. 

The particles from neighboring layers interact with collective potential of this crystal 
and tend to dispose in its minima. In an ideal (translational invariant) case the particles 
have to be found under centers of the neighboring layer exactly. The particles from a next 
layer should be found just under positions of the particles from the first one. However, these 
particles repulse too. So, the energy minimum, should be a nontrivial compromise between 
long-range interaction of the particles and the energy of bonding inside the layers. 

Numerical simulation shows that the symmetric position of the particles over the centers 
of the cells of the Wigner crystal is unstable. They are shifted (in a projection onto layer 
plane) in a direction of one of the neighboring particles. In particular, it breaks the geometric 
frustration in the collective position of a hexagonal 2D lattice from one layer under another. 

It looks like as follows. Moving in a collective field of the neighboring layers the carriers 
attract to the minima of the potential. These minima play a role of effective "positively 
charged particles" . The projections of real particles from the neighboring layers (attracted 
to these minima) on the xy-plane can be treated as "images" of these effective particles. 
This is essentially 3D picture, but in some sense it is quite close to a stripe formation in 
cuprates and nikelatesuu. It is direct analogy to previously studied process of screening in a 
system containing two kinds of the particles that bind into "pairs" . These pairs are dipolarly 
charged and form, in their turn, the chains of the dipoles&l. 

Let us note once more that the "pairs" mentioned above are formed by particles from 
different layers. So, this pairing is essentially 3D phenomenon. Moreover, we found that the 
variation of the density leads to other nontrivial structures also. In particular, at low density 



the system produces an unusual "droplet phase" . This phase consists of the charged droplets 
with an internal structure. One can believe that the specific dipolar chain structure or struc- 
turalised droplet phase can attitude to actively studied now mechanism of superconductivity 
of novel superconductors. 

II. NUMERICAL SIMULATION 

To restore the structures appearing in the system of moving charges placed into 3D 
compensating background with periodically modulated potential along one of the dimensions 
we apply the Brownian dynamics (BD) technique. This technique is widely applied last 
yearsHliil (in particular, by the authors of present articleOo) to simulate a behavior of 
different systems. It is useful in the cases when the discrete description in terms of the 
moving particles is more natural than the continuous kinetic approach. 

The technique is based on the solution of the system of the dynamic equations for the 
particles with discrete set of the coordinates R 3 = {xj,yj,Zj} where 1 < i < iV.and the 
vectors Hjk = Rj — R& = {xj — Xj, yj — yj, Zj — Zj} connecting different particles. The choice 
of appropriate boundary conditions depends on the specific problem. 

Mutual relation between this and kinetic approaches is obtained due to following two 
features of the description. First of all, the equations in such an approach contain some 
random noise sourcesE^Tfcl Together with a relaxation terms in the equations this repro- 
duces an effect of final (nonzero) temperature (in according to known fluctuation-dissipative 
theorem). Besides, if it is necessary, the continual densities of all values of the problem are 
restored "a posteriory" by means of a summation over realizations and averaging during 
sufficiently long time. This time is also determined from the numerical simulation for a 
specific problemists. 

In particular the density of particle distribution in real space Qj(R.) =<< Z)j^(R — 
Rj) >>| to is restored by a summation J2i over particle coordinates and averaging << 
... >>|t over characteristic time t . This density is taken to calculate the correlation func- 



tions of the problem (for example: two-point correlation function G(R, R') = (^(R)^(R')) 
, where (...) denotes an averaging over ensemble of particles). The set of the correlation 
functions gives complete information about the thermodynamic properties of the system. 

More detail description of the approach, as well as some examples of its utilization, can be 
found in previous publications by different authorg9ii2l. In present work , we limit ourselves 
by relatively short positing of the problem and discussion of numerical results. 

The set of BD equations can be written in a form 
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Here 7 is relaxation constant. To model a thermal bath we apply the Gaussian random force 
SF{t), 

(<5F(R;t)>=0; 
(SF(R; t), <5F(R'; t')) = 2 7 T5(R - R')6(t - t') 

Periodic one-particle potential is chosen in following form 

V(Rj) = VoCos^nkRj/Lz) 

The system is supposed to be continued periodically. Here L z /k is a period of the 
potential along the z-axis and k is an integer number. Along two other axes the system 
is quasi-infinite. The boundary conditions along all the axes are as follows. Every particle 
interacts with the particles inside the calculation volume Q = L x ■ L y ■ L z as well as with all 
their "images" obtained as a result of the translation of particle coordinates across (beyond) 
nearest boundaryl3lJ . So, the many-body potential U(Rjk) in the equations contains the 
vectors Xjk = Xjk ± L x , y^ = y^ ± L y , Zjk = Zjk ± L z with all possible permutations. The 
sign "+" or "-" is determined by a direction to the nearest boundary plane along given axis. 
Besides, if moving particle leaves the volume Q = L x ■ L y ■ L z along one of the axes it is 
returned to the volume. At absence of external field it can by done or by means of mirror 
reflection, or by cyclic shift. In the second case corresponding projection of the velocity 



conserves and the coordinate projection is shifted as follows: Xj — > Xj ± L x ,yj — > yj ± L y , 
Zj -► Zj ± zJH. 

This approach is quite suitable for the simulation of translation invariant system 
of charges in solid state. Used here many body potential U(Rjk) = U cou i om b(Tljk) • 
UscreeniRjk) contains from long-range coulomb interaction U cou [ om b(Tljk) ~ 1/ | Rjfc | and 
cut-off screening impact U screen (TLjk)~ exp(— | Hjk | /r ) from crystal lattice and other 
subsystems of the problem. 

Numerical results are summarized in the Figs. [lf|| For definiteness all instant configu- 
rations of the particles here were obtained in the frames of uniform kinetic scenario. The 
random distributions of the charges over system have been taken as the initial conditions. 
From physical point of view appearance of such carriers in the system can be treated as a 
result of a doping. In its turn, it corresponds to standard initial conditions for a search of 
an energy minimum by relaxation techniqueo and annealing techniquefl'Ej. 

Besides, equal number of the particles A=256 has been used for all the cases. A number 
of the variables of 3D problem (the coordinates and velocities) corresponding to this number 
of particles is equal to 2dN = 1536 ~ 1.5 • 10 3 . In its turn, it corresponds to approximately 
2.5 • 10 6 of their mutual combinations that should be calculated at every step. At given 
number of the particles A=256 the charge density is varied by change of the cross-section 
L x ■ L y of the calculation volume Q = L x - L y - L z (at periodic boundary conditions described 
above). To control the results were reproduced at some densities for other numbers of 
particles (N = 128 and iV = 512 , with respectively chosen cross-section of the box L x -L y ). 

III. RESULTS AND DISCUSSION 

At initial stage of the evolution the particles spontaneously group to the vicinities of 
the planes corresponding to the minima of periodic potential V(Rj-) = Vo cos(2tt kHj/L z ). 
These minima below are named as the "layers". For definiteness, we are denote the particles 
as belonging to a given layer in the case if the ^-coordinate is away from the layer not more 
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than 0.1 of distance between the layers \ z — Zk\< 0.1- | Zk+i — Zk |. 

Typical projection on xy-plane of an instant configuration of two arbitrary neighboring 
layers at early evolution stage is shown in the Fig.|l] . The particles from different layers 
are shown by two kinds of the circles (black and white respectively). All other particles are 
not shown in this figure. Typical for this stage tendency to a "pairing" (in projection onto 
xy-pl&ne) of the particles from neighboring layers is seen directly. 

Below this tendency will be characterized numerically. Let us note that such a state 
conserves for some of intermediate densities only. In particular, the configuration shown in 
Fig. [I] is interesting because it takes place in system with relatively high density (L x = L y = 
9 at iV = 256) at initial kinetic stage. Another, tetragonal, structure is found to be natural 
for this system in a final equilibrium state. 

Instant particle configurations at final (but not equilibrium) stages of structure forming 
is shown in Fig. |2| for three specific densities. As before, the only particle coordinates for two 
consistent layers are shown by means of different circles in projection onto xy-plane. The 
case (a) corresponds to already mentioned high density (L x = L y = 9). The configuration 
(b) conserves at long time for the density close to some critical one (L x = L y = 12). And at 
small density (L x = L y = 50) a kind of clasterisation into "droplets" happens in the system. 
These droplets contain of scraps of charge chains and, therefore, possess some kind of fine 
structure. The droplets from different layers repulse mutually. The state obtained looks like 
analogous distributions obtained at concentration stratification (or phase separation) of the 
particles of different kindB. For briefness it could be called as "phase separation" too. 

One can use a time dependence of minimal distance between particles b =< min{| Hj^ | 
} > as a numerical characteristics of different scenarios at different densities g. This value 
corresponds to mean distance between nearest neighbors. Different mutual relations and 
sign of inequality between values of b inside the layers and for the pairs of neighboring layers 
reflects a difference between visually observed structures. 

Four qualitatively different scenarios of evolution of the distance between nearest particles 
are shown in Fig|| It is done for following cases: (a) at high density (L x = L y = 9); (b) 



at some density close, but slightly higher than a critical one (L x = L y = 10.5); (c) (b) at 
a density close, but lower than a critical one (L x = L y = 11); and finally (d) at very low 
density [L x — L y — 70). The averaged projections of the distances inside the layers and 
between them are shown by black and white circles respectively. 

At fixed distance between the layers the average interaction between the layers is pro- 
portional to the density. When the density is high enough this interaction is quite strong to 
order particles in 3D structure. The Fig.f| presents typical tetragonal crystal structure that 
is obtained at high density. The particular case of (L x = L y = 9) is shown in isometric pro- 
jection. Long-wave z-displacements of the particles inside each of layer are directly visible. 
It is caused by strong interaction along z-axis which does not allow to create an ideal 3D 
tetragonal lattice. 

When the density goes to some critical one (10.5 < L x = L y < 11) the force bonding 
the particles in the layers prevail the ^-component of the interaction between layers. The 
particles are locked strongly in the layers and ordered hexagonally. Naturally, the mean 
distance between nearest neighbors stabilizes with a time. However, the distance between 
the projections of particles from different layers does not tend to zero now (as it takes place 
for 3D tetragonal lattice). This fact is reflected by difference in behavior of the lower lines 
in the Figs.||(a-c). 

Finally, when the density is extremely low the repulsion between the particles from 
different layers is notable for sufficiently big groups of the particles. Phase separation occurs 
now. The particles combine into clusters (or droplets). The projections of the clusters from 
different layers mutually fill in free spaces for each other. A tendency to form stripes or 
chains of charge along one of spontaneously chosen direction is observed in the system. 

All these modifications of the structure are reflected by physically measurable corre- 
lation function G(R, R') = (q(R)q(R')) shown in the Fig.|(a-c). The value G{R, R) = 
(g(R)g(R')) is calculated here for one representative layer and shown for following three 
typical cases: (a) well pronounced tetragonal lattice (at L x = L y = 7); (b) for hexagonal or- 
dering of the particles inside the layers [L x = L y = 10.5); (c) for droplets phase appearing at 
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phase separation into charge rich and charge poor domains (at small density L x = L y = 50). 
The well defined filament structure caused by the scraps of charge chains are visible. In 
all pictures the averaged minimal distance between particles inside the layer for the given 
L x = L y are marked by the arrows. 

A dependence of equilibrium averaged minimal distance between the particles b =< 
min{| Rjfc |} > from density (shown as a function from box size L x = L y at fixed number 
of particles N = 256) is summarized in Fig. §. As before, two kinds of the circles denotes 
the values of b corresponding to the particles from the same or different layers. Sloping 
dash-doted line gives a slope of the b(L x ) function for ideal tetragonal lattice. 

This plot can be treated as a phase diagram of the system. This diagram is done in 
terms of the mean density g~l/L x L y at fixed L z . 

Vertical dashed line (at L x m 11) denotes a transition between two different structures. 
It can be proven that for this case the relation between two distances (on upper and lower 
curves) is equal exactly to the a/3. This relation corresponds to ideal 3D hexagonal ordering. 

Dotted vertical line corresponds to a transition from a charge order like Wigner crystals 
to the phase separation structure which contains charge-rich droplets with an internal fine 
structure. It should be noted that the curves b(L x ) do not reach horizontal asymptotes. 

It means that the internal density of the charges in the droplets depends on the value 
of averaged density in the whole system. In spite of some analogy with a picture of liquid 
droplets this fact is caused by absent of direct attraction among particles. The particles can 
not collect together in compact droplets. The only effective attraction exists here which is 
caused by a pairing of those particles from different layers which are close to the droplet 
boundaries. As a result, the size of the droplets is always comparable with a distance between 
them. From physical point of view the dependence of the internal density of charge in the 
droplets from the doping level is one of the mostly important feature of this new droplet 
phase. 
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FIG. 1. Typical instant intermediate particle configuration for two neighbouring layers, pro- 
jected into xy-plane. The particles from different layers are shown by different circles (black 
and white respectively). High density case corresponding to tetragonal equilibrium structure 
(L x = L y = 9 at N = 256) is presented. 
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FIG. 2. Instant configurations of the particles at final stage of evolution are shown for the 



three characteristic densities: (a) at L x = L y = 9 ; (b) at L x = L y = 12 ; (c) at L x = L % 



50 



(TV = 256). The projections of the particles from different layers are shown by black and white 
circles respectively. 
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FIG. 3. Four qualitatively different evolution scenaria of averaged minimal distance between 
the projections of particles from the same layer (black circles) and nearest layers (white circles): 
(a) at L x = L y = 9 ; (b) L x = L y = 10.5;(c) at L x = L y = 11 ; (d) at L x = L y = 70 (N = 256). 
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FIG. 4. Typical tetragonal 3D structure appearing at high density ( L x = L y = 9 ). 
Long-range waves of z-displacements of particles caused by a strong inter-layer interaction are 
seen directly. 
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FIG. 5. Two-point correlation function G(R, R') = (q(R)q(R')) is presented for following 

three typical cases: (a) well pronounced tetragonal lattice (at L x = L y = 7); (b) for hexagonal 
ordering of the particles inside the layers {L x = L y = 10.5); (c) for droplet phase with an internal 
stripe structure appearing at phase separation (at small density L x = L y = 50). The length of 
arrows indicate the value of the averaged minimal distance between particles for given phase. 
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FIG. 6. A dependence of equilibrium averaged minimal distance between the particles 
b =< min{| Rj^ |} > from density (shown as a function from box size L x = L y at fixed number of 
particles N = 256). Two kinds of the circles denotes the values of b corresponding to the particles 
from the same layers (black circles) or different layers (white circles). Sloping dash-doted line gives 
a slope of the b(L x ) function for an ideal tetragonal lattice. Vertical dashed line (at L x ~ 11) 
denotes a transition from tetragonal to hexagonal ordering. Dotted vertical line corresponds to a 
transition from the ordering of separate charges to the concentration stratification. 
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